#ifndef _MIKE_INCLUDE_
#define _MIKE_INCLUDE_

/**
 * @file suborder_relations.hpp
 * @author Mike Jacobson
 * @version $Id$
 */

#include <fstream>
#include <NTL/ZZ.h>
#include <ANTL/quadratic/qo_relation.h>
#include <ANTL/quadratic/qo_partial_relation.h>
#include <ANTL/quadratic/qo_relation_generator.h>

#include<sstream>
#include<string>
#include <stdio.h>
#include <stdlib.h>

#include <givaro/givinit.h>
#include <givaro/givintfactor.h>
#include <givaro/givtimer.h>
#include <givaro/givintprime.h>

extern "C" {
  //#include <stdio.h>
  //#include <stdlib.h>
  //#include <math.h>
  //#include <time.h>
  //#include <string.h>
  //#include <gmp.h>
  //#include <unistd.h>
  //#include <sys/types.h>
  //#include "mpzqform.h"
  //#include "mpzsmoothscan.h"
  #include "smoothrelation.h"
  }

using namespace ANTL;
using namespace std;

void find_ideal_prime_to(const ZZ & a, const ZZ & b, const ZZ & D, const ZZ & f, ZZ & A, ZZ & B) {
    // find equivalent ideal prime to f
    if (GCD(a,f) > 1) {
      ZZ c = (b^2 - D) / (a << 2);
      if (GCD(c,f) > 1) {
	A = a+b+c;
	B = -b - (a << 1);
      }
      else {
	A = c;
	B = -b;
      }
    }
    else {
      A = a;
      B = b;
    }
}


void
go_to_sub_order(const ZZ & a, const ZZ & b, const ZZ & D, const ZZ & f, ZZ & A, ZZ & B)
{
  find_ideal_prime_to(a,b,D,f,A,B);
  mul(B,B,f);
  rem(B,B,A << 1);
}


/*
void
go_to_super_order(const ZZ & a, const ZZ & b, const ZZ & D, const ZZ & f, ZZ & A, ZZ & B)
{
  find_ideal_prime_to(a,b,D*f*f,f,A,B);

  ZZ d,mu,lambda;
  XGCD(d,mu,lambda,f,A);

  mul(B,B,mu);
  if (IsOdd(D)) {
    mul(lambda,lambda,A);
    add(B,B,lambda);
  }
  rem(B,B,(A << 1));
}
*/


void
get_FB_in_Dsub(const ZZ & Dsub, const ZZ & f, factor_base<ZZ> FB_D, qi_class<ZZ> *FB, int sizeFB) \
{
  // set to sub-order
  qi_class<ZZ>::set_current_order(Dsub);

  for (int i=0; i<sizeFB; ++i) {
    ZZ a,b;
    a = FB_D.get_a(i);
    b = FB_D.get_b(i);

    // convert to suborder
    ZZ A,B;
    go_to_sub_order(a,b,Dsub/(f*f),f,A,B);

    // attempt to assign to a prime ideal
    qi_class<ZZ> P;
    if (P.assign(A,B) && f % a != 0)
      FB[i].assign(P);
    else
      FB[i].assign_one();
  }
}


/*
void
get_FB_in_D(const ZZ & D, const ZZ & f, factor_base<ZZ> FB_Dsub, qi_class<ZZ> *FB, int sizeFB) \
{
  // set to sub-order
  qi_class<ZZ>::set_current_order(D);

  for (int i=0; i<sizeFB; ++i) {
    ZZ a,b;
    a = FB_Dsub.get_a(i);
    b = FB_Dsub.get_b(i);

    // convert to suborder
    ZZ A,B;
    go_to_super_order(a,b,D,f,A,B);

    // attempt to assign to a prime ideal
    qi_class<ZZ> P;
    if (P.assign(A,B))
      FB[i].assign(P);
    else
      FB[i].assign_one();
  }
}
*/

// returns the list where list[i]=1 iff prime at index p splits for D1
void 
make_split_list ( vec_long& list_split , const vec_long& p_list , const ZZ& D1 , const ZZ& D2)
{
  list_split.SetLength(0);
  list_split.SetLength(p_list.length());
  for ( int i = 0 ; i < p_list.length() ; i++ )
  {
    if ((Kronecker(D1,p_list[i]) == 1) && (Kronecker(D2,p_list[i]) == 1))
      list_split[i] = 1;
    else
    {
      list_split[i] = 0;
     // cout << p_list[i] << " is non split " << endl;
    }
  }
}

//
// tests whether the vector e is a relation with factor base FB in O(D)
//
bool
is_relation_in_D(const ZZ & D, qi_class<ZZ> *FB, int sizeFB, vec_long & e) {
  // set O(D)
  qi_class<ZZ>::set_current_order(D);

  // compute power-product of prime ideals in O(D2)
  qi_class<ZZ> A;
  A.assign_one();
  for (int i=0; i<sizeFB; ++i) {
    for (int j=0; j<labs(e[i]); ++j) {
	if (e[i] > 0)
          nucomp(A,A,FB[i]);
	else if (e[i] < 0)
	  nucomp(A,A,-FB[i]);
    }
  }
  //  cout << e << endl;
  //  cout << A << endl;

  return A.is_one();
}



//
// Computes the "cardinality" of the relation e in the sense of Gaetan and Drew's paper
//
int
relation_cardinality(const ZZ & D, qi_class<ZZ> *FB, int sizeFB, vec_long & e)
{
  // set to O(D)
  qi_class<ZZ>::set_current_order(D);

  // find indices of, and number of, non-zero entries in e
  int non_zero = 0;
  int indices[sizeFB];
  for (int i=0; i<sizeFB; ++i) {
    if (e[i]) {
      indices[non_zero] = i;
      ++non_zero;
    }
  }
  //  cout << FB[indices[non_zero-1]].get_a() << endl;


  /*
  cout << "non zero indices " ;
  for ( int i = 0 ; i < non_zero ; i++ )
    cout << FB[indices[i]].get_a() << "^" << e[indices[i]] << " ";

  cout << endl;
  */


  vec_long v = e;
  int count = 0;
  int maxidx = 1 << non_zero;

  /*
  cout << "\nnon_zero = " << non_zero << endl;
  cout << "maxidx = " << maxidx << endl;
  */

  qi_class<ZZ> A;
  A.assign_one();
  for (int i=0; i<sizeFB; ++i) {
    for (int j=0; j<labs(e[i]); ++j) {
	if (e[i] > 0)
          nucomp(A,A,FB[i]);
	else if (e[i] < 0)
	  nucomp(A,A,-FB[i]);
    }
  }

  if (A.is_one())
  {
    cout << v << endl;
    ++count;
    //cout << "count = " << count << endl;
  }


  for (int i=1; i<maxidx; ++i) { // Mike says i=0

    //cout << v << endl;


    // update sign array (using Gray code)
    /*
    long j = i;
    int nu_2 = 0;
    while (j>0 && (j & 1) == 0) {
      j >>= 1;
      ++nu_2;
    }
    */

    /*
    int nu_2;
    {
      int two_i = (i << 1);
      int pow_2 = 1;
      int nu    = 0;
      int j;


      while ( !(two_i & 1) ) {
        two_i >>= 1;
        pow_2 <<= 1;
        nu ++ ;
      }

      j = i >> (nu-1);
      if ( j & 1 ) {
        j +=  1;
        j >>= 1;
      }
      else j >>= 1;

      nu_2 = nu-1;
    }
    */

    long j = i;
    int nu_2 = 0;
    while ((j & 1) == 0) {
      j >>= 1;
      ++nu_2;
    }


    //cout << " updating exponent " << indices[nu_2] << "|";

    qi_class<ZZ> B = FB[indices[nu_2]];
    for (int k=1; k<labs(v[indices[nu_2]]); ++k)
      nucomp(B,B,FB[indices[nu_2]]);
    nudupl(B,B);

    if (v[indices[nu_2]] < 0)
      nucomp(A,A,B);
    else
      nucomp(A,A,-B);

    v[indices[nu_2]] *= -1;

    if (A.is_one())
    {
      cout << v << endl;
      ++count;
      //cout << "count = " << count << endl;
    }

  }

  //cout << endl;

  return count;
}


int
relation_cardinality_biasse(const ZZ & D, qi_class<ZZ> *FB, int sizeFB, vec_long & e)
{
  // set to O(D)
  qi_class<ZZ>::set_current_order(D);

  // find indices of, and number of, non-zero entries in e
  int non_zero = 0;
  int indices[sizeFB];
  for (int i=0; i<sizeFB; ++i) {
    if (e[i]) {
      indices[non_zero] = i;
      ++non_zero;
    }
  }

  //  cout << FB[indices[non_zero-1]].get_a() << endl;


  /*
  cout << "non zero indices " ;
  for ( int i = 0 ; i < non_zero ; i++ )
    cout << FB[indices[i]].get_a() << "^" << e[indices[i]] << " ";

  cout << endl;
  */


  vec_long v = e;
  int count = 0;
  int maxidx = 1 << non_zero;

  for ( int i = 0 ; i < maxidx ; i++ )
  {
    qi_class<ZZ> A;
    A.assign_one();

    for ( int k = 0 ; k <  non_zero ; k++ )
    {
      for ( int j = 0 ;  j<labs(e[indices[k]]) ; j++ )
      {
        if ( bit(i , k ) )
          nucomp(A,A,FB[indices[k]]);
        else
          nucomp(A,A,-FB[indices[k]]);
      }
    }

    if (A.is_one())
      ++count;

  }

  //cout << endl;

  return count;
}

int
find_relation_Drew ( const ZZ & D1, const ZZ & D2, const ZZ & v, vec_long & e, vec_long & p_list , int FB_size_in , int FB_smooth_bd)
{
  qi_class<ZZ>::set_current_order(D1);


  qi_class<ZZ>::set_current_order(D1);
  qo_relation_generator<ZZ> GEN;
  GEN.use_large_primes(0);
  GEN.set_SieveRadius(100000);
  GEN.set_P_ONCE(15);
  //GEN.init(D1,to_ZZ(0),SIQS,FB_size_in);
  
  GEN.init(D1,to_ZZ(0),v,SIQS,FB_size_in,0 , false);
  int sizeFB = GEN.FB.get_sizefb();


  qi_class<ZZ> FB_D1[sizeFB];
  p_list.SetLength(sizeFB);
  e.SetLength(sizeFB);

  for (int i=0; i<sizeFB; ++i) {
    ZZ a,b;
    a = GEN.FB.get_a(i);
    b = GEN.FB.get_b(i);
    FB_D1[i].assign(a,b);
    p_list[i] = GEN.FB.get_a(i);
  }

  qi_class<ZZ>::set_current_order(D2);
  qi_class<ZZ> FB_D2[sizeFB];
  for (int i=0; i<sizeFB; ++i)
    FB_D2[i].assign_prime(to<ZZ>(GEN.FB.get_a(i)));


  int count_D1,count_D2;
  int seed_drew = 1;
  stringstream str_D , str_bd , str0 , str1;
  char** input_drew = new char*[5];

  str_D  << D1 ;
  str_bd << FB_smooth_bd ;
  str0   << 0;
  str1   << 1;

  input_drew[1] = new char[500];
  input_drew[2] = new char[100];
  input_drew[3] = new char[100];
  input_drew[4] = new char[100];

  strcpy ( input_drew[1] , str_D.str().c_str() );
  strcpy ( input_drew[2] , str_bd.str().c_str() );
  strcpy ( input_drew[3] , str0.str().c_str() );
  strcpy ( input_drew[4] , str1.str().c_str() );
  //cout << "p_list = " << endl << p_list << endl;

  vec_long list_split;
  make_split_list ( list_split , p_list , D1 , D2 );

  do
  {
    int m = test_drew ( 5 , input_drew , (unsigned long)seed_drew);

    //for ( int i = 0 ; i < m ; i++ ) if ( h[i] ) printf ("%ld^%s%ld ", ell[i], (hs[i]<0?"-":""), h[i]);

    for ( int i = 0 ; i < e.length() ; i++ )
      e[i] = 0;


    int all_split = 1;
    seed_drew ++;
    for ( int i = 0 ; i < m ; i++ )
    {
      if ( h[i] )
      {
        int found_p = 0;
        for ( int k = 0 ; k < p_list.length() ; k++ )
        {
          if ( p_list[k] == ell[i] )
          {
            found_p = 1;
            e[k]    = h[i];
          }
        }
        if (!found_p)
          all_split = 0;
      }
    }

    
    for ( int i = 0 ; i < e.length() ; i++ )
    {
      if ( e[i] )
      {
        cout << p_list[i] << "^" << e[i] << " ";
      }
    }
    cout << endl;
    

    qi_class<ZZ>::set_current_order(D1);

    /*
    qi_class<ZZ> B = (qi_class<ZZ>) GEN.FB.product_over_FB(e);
    cout << "Power product:  " << B << endl;
    if (!B.is_one()) {
      cerr << "ERROR:  relation in valid!" << endl;
    }
    */

    // compute relation cardinalities over Dsub and D
   
    // Checks if the relation only involves split primes. 
    // If not, it can't be evaluated in End(E) by Drew's 
    // code.
    
    for ( int i = 0 ; i < e.length() ; i++ )
    {
      if ( e[i] )
      {
        if ( !list_split[i] )
          all_split = 0;
      }
    }


    if ( all_split ) 
    {
      cout << "before count" << endl;
      cout << "D1 = " << D1 << endl;
      cout << "D2 = " << D2 << endl;
      cout << "e = " << e << endl;
      cout << "size_FB = " << sizeFB << endl;

      for ( int k = 0 ; k < sizeFB ; k++ )
        cout << FB_D2[k] << " ";
      cout << endl;

      count_D1 = relation_cardinality(D1,FB_D1,sizeFB,e);
      int count_biasse = relation_cardinality_biasse(D1,FB_D1,sizeFB,e);
      count_D2 = relation_cardinality(D2,FB_D2,sizeFB,e);
      if ( count_D1 == 0 )
      {
        cout << "Relation valid for Drew but rejected by relation_cardinality() !" << endl;
        cout << "count_biasse = " << count_biasse << endl;
      }
      else
      {
        cout << "count_biasse = " << count_biasse << endl;
        cout << "count_D1     = " << count_D1 << endl;

        if ( count_D1 != count_biasse )
        {
          cout << "Mike's count gives something else" << endl;
          count_D1 = count_biasse;
        }
      }
    }
    else
    {
      count_D1 = 0;
      cout << "relation with a non split prime" << endl;
    }


  }while(count_D1 <= count_D2 && count_D1 == 0);

  return count_D1;

}





//
// This is the sieving version of the FindRelation algorithm in Gaetan and Drew's paper.
//

int
find_relation_sieve(const ZZ & D1, const ZZ & D2, const ZZ & v, vec_long & e, vec_long & p_list , int FB_size_in)
{
  // set relation generation parameters
  qi_class<ZZ>::set_current_order(D1);
  qo_relation_generator<ZZ> GEN;
  GEN.use_large_primes(0);
  //  GEN.set_SieveRadius(250000);
  // 	GEN.set_P_ONCE(28);
  GEN.set_ToleranceValue(6.0);
  GEN.use_bernstein(500);
  //  GEN.init(D1,to_ZZ(0),MPQS,340);
  GEN.set_SieveRadius(100000);
  GEN.set_P_ONCE(15);
  //GEN.init(D1,to_ZZ(0),SIQS,FB_size_in);

  GEN.init(D1,to_ZZ(0),v,SIQS,FB_size_in,0 , true);
  int sizeFB = GEN.FB.get_sizefb();

  //cout << "Mike sieve code : sizeFB = " << sizeFB << endl;

  // compute factor base ideals in O(D1)
  qi_class<ZZ>::set_current_order(D1);
  qi_class<ZZ> FB_D1[sizeFB];
  p_list.SetLength(sizeFB);
  for (int i=0; i<sizeFB; ++i) {
    ZZ a,b;
    a = GEN.FB.get_a(i);
    b = GEN.FB.get_b(i);
    FB_D1[i].assign(a,b);
    p_list[i] = GEN.FB.get_a(i);
  }

  //cout << "sieve Mike, max p = " << p_list[sizeFB-1] << endl;

  //cout << endl << p_list << endl;

  // compute factor base ideals in D2
  qi_class<ZZ>::set_current_order(D2);
  qi_class<ZZ> FB_D2[sizeFB];
  for (int i=0; i<sizeFB; ++i)
    FB_D2[i].assign_prime(to<ZZ>(GEN.FB.get_a(i)));

    

  // find relations until a suitable one is found
  int count_D1, count_D2;
  vec_long list_split;
  make_split_list ( list_split , p_list , D1 , D2);

  do {
    // set to super order
    qi_class<ZZ>::set_current_order(D1);

    qo_relation<ZZ> rel;
    long j1,j2;
    GEN.get_relation(rel,j1,j2);
    rel.get_vector(e);

    int all_split = 1;
    for ( int i = 0 ; i < e.length() ; i++ )
    {
      if ( e[i] )
      {
        if ( !list_split[i] )
          all_split = 0;
      }
    }



    if ( all_split )
    {
      //cout << "\nRelation in O(D1):  " << e << endl;
      // make sure it is a relation!
      qi_class<ZZ> B = (qi_class<ZZ>) GEN.FB.product_over_FB(e);
      //cout << "Power product:  " << B << endl;
      if (!B.is_one()) {
        cerr << "ERROR:  relation in valid!" << endl;
      }

      // compute relation cardinalities over Dsub and D
      count_D1 = relation_cardinality(D1,FB_D1,sizeFB,e);
      count_D2 = relation_cardinality(D2,FB_D2,sizeFB,e);
    }
    else
    {
      count_D1 = 0;
      cout << "we have a relation with non-split primes" << endl;
    }
    //cout << "count_D1 = " << count_D1 << ", count_D2 = " << count_D2 << endl;

  } while(count_D1 <= count_D2 && count_D1 == 0);

  return count_D1;
}



//
// finds a relation that holds in O(D) but not O(Dsub), where O(Dsub) \in O(D)
//
void
find_relation_sieve_nocount(const ZZ & Dsub, const ZZ & D, vec_long & e, vec_long & p_list) 
{
  // compute conductor of O(Dsub) relatinve to O(D)
  ZZ f,q,r;
  DivRem(q,r,Dsub,D);
  f = SqrRoot(q);
  if (!IsZero(r) || f*f != q) {
    cerr << "ERROR:  Must have Dsub/D a square" << endl;
    exit(1);
  }

  // set relation generation parameters
  qi_class<ZZ>::set_current_order(D);
  qo_relation_generator<ZZ> GEN;
  //  GEN.set_SieveRadius(250000);
  // 	GEN.set_P_ONCE(28);
  //  GEN.set_ToleranceValue(2.3);
  //	GEN.use_bernstein(200);
  //  GEN.init(D1,to_ZZ(0),MPQS,340);
  GEN.set_SieveRadius(130000);
  GEN.set_P_ONCE(9);
  //GEN.init(D,to_ZZ(0),SIQS,500);

  GEN.init(D,to_ZZ(0),to_ZZ(0),SIQS,500,0 , true);
  // compute factor base ideals in O(Dsub)
  int sizeFB = GEN.FB.get_sizefb();
  qi_class<ZZ>::set_current_order(Dsub);
  qi_class<ZZ> FB_Dsub[sizeFB];
  get_FB_in_Dsub(Dsub,f,GEN.FB,FB_Dsub,sizeFB);


  // find a relation in O(D) but not in O(DSub)
  do {
    // set to super order
    qi_class<ZZ>::set_current_order(D);

    qo_relation<ZZ> rel;
    long j1,j2;
    GEN.get_relation(rel,j1,j2);
    rel.get_vector(e);

    cout << "\nRelation in O(D):  " << e << endl;

    // make sure it is a relation!
    qi_class<ZZ> B = (qi_class<ZZ>) GEN.FB.product_over_FB(e);
    cout << "Power product:  " << B << endl;
    if (!B.is_one()) {
      cerr << "ERROR:  relation in valid!" << endl;
    }

  } while (is_relation_in_D(Dsub,FB_Dsub,sizeFB,e));


  // set list of primes for relation
  qi_class<ZZ>::set_current_order(D);
  p_list.SetLength(sizeFB);
  for (int i=0; i<sizeFB; ++i)
    p_list[i] = GEN.FB.get_a(i);
}

int
FB_size_from_FB_bound ( const ZZ& Delta_in , long bound_p )
{
  PrimeSeq primes;        //generates primes
  long p;        
  int  count;
  ZZ   temp , Delta;

  p     = primes.next();
  count = 0;
  Delta = Delta_in;

  while ( p < bound_p )
  {
    p = primes.next();

    if (Kronecker(Delta,p) != -1)
    {
      rem(temp, Delta, to_ZZ(p*p));
      if (temp != 0)
        count ++;
    }
  }

  return count;
}

/* Finds n relations without counting or anything else */
void
find_relation_sieve_biasse(
    const ZZ & D1, 
    const ZZ & D2, 
    const ZZ & v, 
    vec_long & e, 
    vec_long & p_list , 
    int FB_size_in ,
    double Tol ,
    long  sieve_r , 
    int   P_once ,
    int rel_nb)
{
  // set relation generation parameters
  qi_class<ZZ>::set_current_order(D1);
  qo_relation_generator<ZZ> GEN;
  GEN.use_large_primes(0);
  //  GEN.set_SieveRadius(250000);
  // 	GEN.set_P_ONCE(28);
  GEN.set_ToleranceValue(Tol);
  //	GEN.use_bernstein(200);
  //  GEN.init(D1,to_ZZ(0),MPQS,340);
  GEN.set_SieveRadius(sieve_r);
  GEN.set_P_ONCE(P_once); //6
  //GEN.init(D1,to_ZZ(0),SIQS,FB_size_in);

  GEN.init(D1,to_ZZ(0),to_ZZ(0),SIQS,FB_size_in,0 , true);
  int sizeFB = GEN.FB.get_sizefb();

  //cout << "Mike sieve code : sizeFB = " << sizeFB << endl;
  //cout << "FB_size_in was " << FB_size_in << endl;
#if 0
  // compute factor base ideals in O(D1)
  qi_class<ZZ>::set_current_order(D1);
  qi_class<ZZ> FB_D1[sizeFB];
  p_list.SetLength(sizeFB);
  for (int i=0; i<sizeFB; ++i) {
    ZZ a,b;
    a = GEN.FB.get_a(i);
    b = GEN.FB.get_b(i);
    FB_D1[i].assign(a,b);
    p_list[i] = GEN.FB.get_a(i);
  }

  cout << p_list[sizeFB-1];


  // compute factor base ideals in D2
  qi_class<ZZ>::set_current_order(D2);
  qi_class<ZZ> FB_D2[sizeFB];
  for (int i=0; i<sizeFB; ++i)
    FB_D2[i].assign_prime(to<ZZ>(GEN.FB.get_a(i)));
#endif
    

  // find relations until a suitable one is found
  int count_D1, count_D2;

  count_D1 = 0;
  do {
    // set to super order
    qi_class<ZZ>::set_current_order(D1);

    qo_relation<ZZ> rel;
    long j1,j2;
    GEN.get_relation(rel,j1,j2);
    rel.get_vector(e);

    for ( int i = sizeFB-1 ; i >= 0 ; i-- )
    {
      if ( e[i] != 0 )
        cout << GEN.FB.get_a(i) << "^" << labs(e[i]) << " ";
    }

    count_D1 ++;

#if 0
    cout << "\nRelation in O(D1):  " << e << endl;
    // make sure it is a relation!
    qi_class<ZZ> B = (qi_class<ZZ>) GEN.FB.product_over_FB(e);
    cout << "Power product:  " << B << endl;
    if (!B.is_one()) {
      cerr << "ERROR:  relation in valid!" << endl;
    }

    // compute relation cardinalities over Dsub and D
    //count_D1 = relation_cardinality(D1,FB_D1,sizeFB,e);
    //count_D2 = relation_cardinality(D2,FB_D2,sizeFB,e);

    //    cout << "count_D1 = " << count_D1 << ", count_D2 = " << count_D2 << endl;
#endif
  } while(count_D1 < rel_nb);

}

int
relation_cardinality_end_E ( const vec_long& e_vec , const vec_long& p_vec , const ZZ& p , const ZZ& j)
{
  stringstream str_command;
  char* count_char = new char[100];

  str_command << "./zp_count_relation " << p << " " << j << " ";
  for ( int i = 0 ; i < e_vec.length() ; i++ )
    str_command << p_vec[i] << "^" << e_vec[i] << " ";

  str_command << " > output_mod_pol ";
  int temp = system ( str_command.str().c_str() );

  while ( !WIFEXITED(temp) ){}

  ifstream myfile("zp_count_relation.log");
  if ( myfile.is_open())
  {
    int k = 0;
    string line;
    getline (myfile,line);

    while ( line[k] != ' ')
    {
      count_char[k] = line[k];
      k ++;
    }

    count_char[k] = '\0';
    int count     = atoi(count_char);

    system ( "rm zp_count_relation.log");
    return count;
  }
  else
  {
    cout << "error relation cardinality : couldn't open log" << endl;
    return 0;
  }






  /*
  ifstream myfile("output_mod_pol" );
  if ( myfile.is_open())
  {
    int found = 0;
    while ( !found ) 
    {
      string line;
      getline (myfile,line);

      //cout << "we have the line " << line << endl;

      if ( line[0] == 'R' )
      {
        found = 1;
        int k = 15;
        while ( line[k] != ' ' )
        {
          count_char[k-15] = line[k];
          k ++;
        }

        count_char[k-15] = '\0';
      }
    }

    int count = atoi(count_char);
    //system ( "rm output_mod_pol" );

    return count;
  }
  else
  {
    cout << "Error count relation in End(E) : couldn't open the file" << endl;
    return 0;
  }
  */
}

void factor_v ( const ZZ& v , vec_long& vec_e , vec_ZZ& vec_p )
{
  stringstream str_v;
  std::vector<Integer> Lf,  Lo;
  IntFactorDom<> IP;

  Integer m;

  str_v << v;
  m = Integer(str_v.str().c_str());

  IP.set(Lf ,Lo , m , 0);

  vec_p.SetLength(0);
  vec_e.SetLength(0);

  for ( int i = 0 ; i < Lo.size() ; i++ )
  {
    stringstream str_p , str_e;
    ZZ p ;
    long e;

    str_p << Lf[i];
    str_e << Lo[i];

    str_p >> p;
    str_e >> e;

    append ( vec_p , p );
    append ( vec_e , e );
  }
}





#endif
